import numpy as np
import pandas as pd
import math
import statsmodels.api as sm
import matplotlib.pyplot as plt
from sklearn.decomposition import PCA
from sklearn.preprocessing import scale
from sklearn.impute import SimpleImputer
pd.set_option("display.max_columns",500)
pd.options.mode.chained_assignment = None  # default='warn'

#所有的表都放在当前文件夹的子文件夹res中
# 原始数据表是：2020data01.xlsx
# 最终主成分分析后得到的表是：fiannly.xlsx
# 其他表为中间结果表，用于观察数据
#批量注释键：Ctrl+/,按这个组合键也可以对已经注释的内容解注释

def newHpFilter(path):
    initExcel = pd.read_excel(path, sheet_name='initial', header=None)
    # 数值部分：[5:,1:]
    initData = initExcel.loc[5:, 1:]  # initData的行标是和initExcel一样的，但是列标重新从1到76编号
    # print(initData)
    hpT = pd.DataFrame()  # 用来装载hp滤波后的dataframe
    subMat=pd.DataFrame() #装cycle的
    notNa = ~initData.isna()
    for col in range(1, initData.shape[1] + 1):
        cycle, trend = sm.tsa.filters.hpfilter(initData.loc[notNa.loc[:, col], col],lamb=1600)
        subMat=pd.concat([subMat,cycle],axis=1,ignore_index=True)
        hpT = pd.concat([hpT, trend], axis=1,ignore_index=True)
    #return subMat #行标：5:76 列标：1：2159
    stdMat=subMat.std(ddof=0,axis=0)
    newSubMat=subMat.copy(deep=True)
    for ncol in range(0,subMat.shape[1]):
        for nrow in range(0,subMat.shape[0]):
            if subMat.iloc[nrow,ncol]>1.5*np.abs(stdMat.iloc[ncol]):
                newSubMat.iloc[nrow, ncol] = 1.5 * np.abs(stdMat.iloc[ncol])
            elif subMat.iloc[nrow,ncol]<-1.5*np.abs(stdMat.iloc[1]):
                newSubMat.iloc[nrow, ncol] = -1.5 * np.abs(stdMat.iloc[ncol])


    # newSubMat.to_excel(excel_writer='D:\\subMat.xlsx')
    # hpT.to_excel(excel_writer='D:\\hpT.xlsx')
    resMat=newSubMat+hpT #resMat 任然是index:[5,76]  columns:[0,2158] notNa:index:[5:76],columns:[1:2159]

    # print(resMat)
    resMat.columns=notNa.columns
    fAns=pd.DataFrame()
    # print(notNa)
    # print(resMat)
    for rCol in range(1, initData.shape[1] + 1):
        Cycle1, Trend1 = sm.tsa.filters.hpfilter(resMat.loc[notNa.loc[:, rCol], rCol], 10)
        fAns=pd.concat([fAns,Trend1],axis=1)
    #fAns.to_excel(excel_writer='D:\\shuju.xlsx',sheet_name='sheet1',header=None,index=False)
    fAns.columns=resMat.columns
    initExcel.loc[5:, 1:] = fAns
    # 交换第0行和第4行位置，让X1,X2,X3...这一行变成第一行
    a, b = initExcel.loc[0, :].copy(deep=True), initExcel.loc[4, :].copy(deep=True)
    initExcel.loc[0, :], initExcel.loc[4, :] = b, a
    # 处理日期格式，不然写入excel后日期就变成了#######
    initExcel.loc[5:, 0] = initExcel.loc[5:, 0].apply(lambda x: str(pd.to_datetime(x).date()))
    # 指定header=None，Index=False，不然会将行、列的索引号也写入excel中
    initExcel.to_excel(excel_writer='.\\res\\in新的python跑出来hp滤波.xlsx', sheet_name='Sheet1', header=None, index=False)
    return initExcel


#cross_cor
# 测试代码：
# x=pd.DataFrame(np.array([3.98,8.23,9.9,3.3,2.95,21.43,17.64,8.222]))
# y=pd.DataFrame(np.array([1.38,2.90,4.34,3.64,5.89,7.56,8.23,9.53]))
# L=6
# ansT=cross_cor(x,y,L)
# print(ansT)
#经过和matlab的多组不同数据结果对比，该函数正确
def cross_cor(x,y,L)->pd.DataFrame:#x,y都是Series
    print('到cross_cor函数了')
    T=y.shape[0]
    p=pd.DataFrame(np.zeros((1,2*L+1)),dtype=float)#1行，2L+1列
    varx=(x-x.mean())*(x-x.mean())
    print('varx is'+str(varx))
    # sigma_x=varx.sum().apply(lambda x:x**0.5)[0]#sum()返回的是一个1*1的df，apply后还是1*1
    sigma_x= np.sqrt(varx.sum()[0])#sum()返回的是一个1*1的df，apply后还是1*1
    vary=(y-y.mean())*(y-y.mean())
    # sigma_y = vary.sum().apply(lambda y: y ** 0.5)[0]
    sigma_y = np.sqrt(vary.sum()[0])
    print(sigma_x)
    print(sigma_y)
    j=0
    for i in range(-L,0):#取不到0
        n1=T+i
        n2=-i
        u=(y.iloc[0:n1,:]-y.mean())
        u=u.values
        v=x.iloc[n2:,:]-x.mean()
        v=v.values
        uv=u*v
        # m1=((y.iloc[0:n1,:]-y.mean())*(x.iloc[n2:,:]-x.mean()))#这样写不行，因为(y.iloc[0:n1,:]-y.mean())和x.iloc[n2:,:]-x.mean()的index不同，*会变成Nan
        s1=uv.sum(axis=0)
        # ss1=s1[0]
        p[j]=pd.DataFrame([s1[0]/(sigma_y*sigma_x)])
        # print(p[j])
        j+=1
    for i in range(0,L+1):
        n1=T-i
        n2=i
        e=(x.iloc[0:n1,:]-x.mean())
        e = e.values
        r=(y.iloc[n2:,:]-y.mean())
        r=r.values
        er=e*r
        s2=er.sum(axis=0)
        # print(s2[0])
        p[j]=pd.DataFrame([s2[0]/(sigma_x*sigma_y)])
        # print(p[j])
        j+=1
    period=pd.DataFrame(np.arange(L,-L-1,-1)).T
    p=pd.concat([period,p],axis=0,ignore_index=True)#纵向拼接
    p=p.T
    p=pd.concat([p,p[1].abs()],axis=1,ignore_index=True)#增加一列abs值，方便查找至最大的那一列，待会不取就完事
    #现在的p第0列是期数，第1列是对应的相关系数，第2列是第1列的绝对值，要找第2列的最大值的索引
    # print(p)
    #新的
    # p=p.fillna(value=np.NINF)
    p = p.values
    idx = p[:, 2].argmax(axis=0)
    # print(idx)
    ansT = p[idx, [1, 0]]
    ansDf = pd.DataFrame(ansT)
    ansDf = ansDf.T
    # print(ansDf)
    return ansDf

#求解相关系数
def getCorre(path,L):
    print('这里是getCoree函数')
    excel=pd.read_excel(path,sheet_name='Sheet1',header=None)#header=None:代表不使用源表的表头
    numData=excel.iloc[5:,1:]
    numData=numData.reset_index(drop=True)#注意，现在的numData行索引从0开始，列索引从1开始，因为没有重置列索引
    print(numData)
    [h,k]=numData.shape
    rho=pd.DataFrame(np.zeros((k,2))*np.nan)
    xnan=~numData.isnull()
    for col in range(1,k+1):#[1,k]，能取到k的
        #nx的index还是原来的，即如果非空的行为12,13...，那么nx的第一行的index就是12
        nx=numData.loc[xnan[col],col].values#iloc是按照索引来的，索引[:,1]取的是第2列而不是第一列，而loc[:,1]就是取得第一列
        ny=numData.loc[xnan[col],1].values#注意这里ny也要取xnan[col]，这样才是对应的
        nx=(nx-np.mean(nx))/np.std(nx,axis=0)
        ny=(ny-np.mean(nx))/np.std(ny,axis=0)
        x=pd.DataFrame(nx)
        y=pd.DataFrame(ny)
        res=cross_cor(x,y,L)
        # print(res)
        res=res.values
        res=pd.DataFrame(res).loc[0]
        rho.loc[col,:]=res
    rho=rho.dropna()
    rho_sorted=rho.sort_values(by=[1,0],axis=0)
    flag=rho_sorted.index
    print(rho_sorted)
    # print(excel)
    # print(rho_sorted)
    # print(flag)
    name=excel.iloc[0:5,1:].T
    # name=name.reset_index(drop=True)//reindex后变成了0开始，到2158了
    print(name)
    name_sorted=name.loc[flag,:]
    ansDf=pd.concat([name_sorted,rho_sorted],axis=1,ignore_index=True)
    ansDf.to_excel(excel_writer='.\\res\\new求解相关系数.xlsx',sheet_name='Sheet1',header=None,index=False)
    return ansDf
    #测试数据
    # path = 'D:\\pyhpfilter.xlsx'
    # L = 6
    # ansDf = getCorre(path, L)
    # print(ansDf)

#获取每一期绝对值大于0.6且最大的那num个
def getMaxAbs(path,num,L):#path是文件路劲，num是取前num个，L是期
    # path='D:\\ProgramData\\MATLABProjects\\'
    table=pd.read_excel(path,sheet_name='Sheet1',header=None)
    table=table.rename(columns={0:'X标号',1:'指标名称',2:'频率',3:'单位',4:'指标ID',5:'相关系数',6:'期'})
    table['sort']=table['相关系数'].abs()
    table=table.sort_values(by=['期','sort'],ascending=[True,False])
    table=table.drop(labels=['sort'],axis=1)
    table.to_excel(excel_writer='myTestBest.xlsx',header=True,sheet_name="Sheet1")

    ans=pd.DataFrame(columns=['X标号','指标名称','频率','单位','指标ID','相关系数','期'])
    for i in range(-L,L+1):
        tmp=table.loc[(table['期']==i) & (table['相关系数'].abs()>=0.6)]
        # tmp=table[table['相关系数']>=0.6 and (table['期']==i)]
        if tmp.shape[0]>num:
            tmpAns=tmp.head(num)
            ans=pd.concat([ans,tmpAns],axis=0,ignore_index=True)
        else:
            ans=pd.concat([ans,tmp],axis=0,ignore_index=True)
    print(ans)
    ans.to_excel(excel_writer='.\\res\\选取相关度最高的n个.xlsx',header=None,sheet_name='Sheet1',index=None)

def myPCA(path1,path2):#path1是hp滤波后的，path2是相关系数表
    excel1=pd.read_excel(io=path1,sheet_name='Sheet1')#以第一行作为表头
    excel2=pd.read_excel(io=path2,sheet_name='Sheet1',header=None)
    excel2=excel2.rename(columns={0:'X标号',1:'频率',2:'单位',3:'指标ID',4:'指标名称',5:'相关系数',6:'期'})
    group=excel2.groupby('期')#分组聚集，size()方法返回每组的行数
    gSize=group.size()#gSize是一个Series
    gNum=gSize.shape[0]#gNum是组数，等于2*L+1
    L=(gNum-1)//2
    pos=0
    writer = pd.ExcelWriter('.\\res\\fiannly.xlsx')
    for t in range(0,gNum):
        #每组数据的行标：[pos,pos+gSize.iloc[t]),每次执行完后pos=pos+gSize.iloc[t]
        Xids=excel2.iloc[pos:pos+gSize.iloc[t],0]
        XidsT=Xids.T#转置一波，方便直接在excel1里面利用列名查找
        goalCols=excel1[XidsT]#这就是我们要找的这一期的那些列
        for row in range(pos,pos+gSize.iloc[t]):
            if excel2.iloc[row,5]<0:#这个变量的相关系数小于0，那么那一列数据去相反数
                Xid=excel2.iloc[row,0]
                goalCols.loc[4:,Xid]=-goalCols.loc[4:,Xid]
        goalData=goalCols.iloc[4:,:]#取数值部分,index和goalCols一样，从4开始，columns变成指定的列
        notnaMat=goalData.isna()
        notnaMat=notnaMat.reset_index(drop=True)

        # print(notnaMat)
        goalDataNp=goalData.values#转出numpy数组
        imp_mean = SimpleImputer(missing_values=np.nan, strategy='mean')#预处理，填充空值为每一列的均值
        imp_mean.fit(goalDataNp)
        goalDataNp=imp_mean.transform(goalDataNp)
        npGoalData=scale(goalDataNp)#标准化，标准化之后就自动根据协方差矩阵进行主成分分析了
        pca = PCA()  # 可以调整主成分个数，n_components = 1
        pca.fit(npGoalData)
        # print(pca.explained_variance_)  # 输出特征根
        # print(pca.explained_variance_ratio_)  # 输出解释方差比
        # print(pca.components_)  # 输出主成分
        pcaComponent=pca.components_#主成分矩阵
        pcaFirstCol=pcaComponent[:,0]#，只取第一列
        pcaFirstRow=pcaFirstCol.T
        ansNp=pcaFirstRow*goalDataNp
        resDf=pd.DataFrame(ansNp)#此时resDf的index已经从0开始了，原来的是从4开始，双射，逆映射
        notnaMat.columns=resDf.columns
        for cool in range(0,notnaMat.shape[1]):#把之前是空值但被填充为均值的那些数全部变为空值
            resDf.loc[notnaMat.loc[:,cool],cool]=np.nan
        print(notnaMat)
        print(resDf)
        #开始处理文本部分为正确格式
        text=goalCols.iloc[0:4,:]#文本部分
        text=text.T#通过转置删除原来的columns，不然拼接不成功
        text=text.reset_index(drop=True)
        text=text.T
        header=goalCols.columns
        header=pd.DataFrame(header).T
        ansText=pd.concat([header,text],ignore_index=True)
        #至此，成功得到了正确格式的文本部分
        #拼接文本和数值部分，两者的列索引是一样的，所以不用改,直接拼接，然后写入文件
        ans=pd.concat([ansText,resDf],axis=0,ignore_index=True)
        gdpColName=excel1.columns.tolist()[1]
        other=excel1.iloc[0:,1]
        gdpColName=pd.DataFrame([gdpColName])
        gdpCol=pd.concat([gdpColName,other],axis=0,ignore_index=True)
        ans=pd.concat([gdpCol,ans],axis=1,ignore_index=True)
        ans.to_excel(writer,sheet_name='Sheet'+str(t-L),header=True,index=None)
        print('pos='+str(pos))
        pos = pos + gSize.iloc[t]
    writer.save()
    writer.close()
    return

if __name__ == '__main__':
    path='.\\res\\2020data01.xlsx'
    newHpFilter(path)
    getCorre('.\\res\\in新的python跑出来hp滤波.xlsx',4)
    getMaxAbs('.\\res\\new求解相关系数.xlsx',num=10,L=4)
    myPCA('.\\res\\in新的python跑出来hp滤波.xlsx','.\\res\\选取相关度最高的n个.xlsx')
    print('跑完了')